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We propose a scheme for constructing classical spin Hamiltonians from Hunds coupled spin-fermion 
models in the limit Jn/t — > oo. The strong coupling between fermions and the core spins requires 
self-consistent calculation of the effective exchange in the model, either in the presence of inhomo- 
geneities or with changing temperature. In this paper we establish the formalism and discuss results 
mainly on the "clean" double exchange model, with self consistently renormalised couplings, and 
compare our results with exact simulations. Our method allows access to system sizes much beyond 
the reach of exact simulations, and we can study transport and optical properties of the model 
without artificial broadening. The method discussed here forms the foundation of our papers Phys. 
Rev. Lett. 91, 246602 (2003), and Phys. Rev. Lett. 92, 126602 (2004). 



I. INTRODUCTION 



The double exchange (DE) model was introduced by 
Zener 1 in 1951 to motivate ferromagnetism in the per- 
ovskite manganites. In contrast to 'Heisenberg like' cou- 
pling between localised spins, the effective interaction 
in 'double exchange' arises from optimisation of carrier 
kinetic energy in the spin background. The intimate 
correlation between spin configuration and electron mo- 
tion had, till recently, restricted the study of the DE 
model to mostly qualitative analysis or mean field the- 
ory. The original proposal of Zener was followed up 2 by 
Anderson and Hasegawa, who clarified the physics of the 
coupled spin-fermion system in a two site model, and 
de Gennes 3 who presented a thermodynamic calculation 
and a phase diagram (incorporating antiferromagnetic 
superexchange) . He produced the first estimate of transi- 
tion temperature (T c ) in the model. The thermodynamic 
transition within double exchange was also studied 4 by 
Kubo and Ohata. This short list essentially exhausts ac- 
tivity on the double exchange problem prior to the 'man- 
ganite renaissance'. 

The discovery of colossal magnetoresistance (CMR) 
and a variety of magnetic phases in the manganites 5 led 
to renewed interest in the DE model. In addition, the 
availability of powerful analytical and numerical tools, 
e.g, dynamical mean field theory (DMFT) and Monte 
Carlo methods provided impetus for studying the DE 
model in detail. In real systems the double exchange in- 
teraction is supplemented by 6 antiferromagnetic (AF) su- 
perexchange, electron-phonon interactions, and disorder, 
and some of these models have been studied within var- 
ious approximations. The primary limitation of current 
methods, as we discuss in detail later, is their inability 
to access transport properties taking spatial fluctuations 
and disorder effects fully into account. In this context 
our method, of constructing an approximate but explicit 
classical spin Hamiltonian, allows a breakthrough. In 
the present paper our detailed results are on the simplest 



case, of the clean DE model. In earlier short publica- 
tions we have presented results on the disordered double 
exchange model 7 , and on magnetic phase competetion 8 . 

Let us define the general model to which our method 
is applicable. H = H e i + Haf, with 



Haf — Js ^ Si.Sj 



(1) 



The U 



-t are nearest neighbour hopping, on a 



square or cubic lattice as relevant, is the on site po- 
tential, uniformly distributed between ±A/2, say, and Js 
is an antiferromagnetic superexchange between the core 
spins. Jh is the 'Hunds' coupling, and we will work in 
the limit Jn/t — » oo. The parameters in the problem are 
A/t, Js/t, and the carrier density n (or chemical poten- 
tial /i). We assume a classical core spin, setting |S;| = 1, 
and absorb the magnitude of the spin in Js- All our en- 
ergy scales, frequency (w) and temperature (T), etc, will 
be measured in units of t. 

For Jh ft — » oo the fermion spin at a site is constrained 
to be parallel to the core spin, gaining energy —Jjj/2, 
while the 'antiparallel' orientation is pushed to +Jh/2. 
Since the hopping term £y itself is spin conserving, the 
motion of the low energy, locally parallel spin, fermions 
is now controlled by nearest neighbour spin orientation. 
The strong magnetic coupling (Jh) generates an effec- 
tive single band 'spinless fermion' problem 9 , with core 
spin orientation dependent hopping amplitudes. We will 
discuss the hopping term further on, for the moment let 
us denote the renormalised (spin orientation dependent) 
hopping amplitude as t, indicative of double-exchange 
physics. 

The t — A — Js problem has a variety of ground states. 
(i) In the absence of Js, both the 'clean' and the dis- 
ordered DE model has a ferromagnetic ground state, at 
all electron density, with T c reducing with increase in 
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A. (ii) The non disordered problem, with J5, leads to 
a variety of phases 10 ' 11 competing with ferromagnetism. 
These are ferromagnetic and A, C, G type AF phases, 
etc. There could also be more exotic 'flux', 'skyrmion' or 
'island' phases in some parts of parameter space. The 
boundaries between these phases are often first order 
so there are regimes of macroscopic phase coexistence. 
The specific set of possible AF phases depends on Jg. 
(Hi) Weak disorder in the t — Js problem 8 ' 12 converts 
the regions of macroscopic phase separation into meso- 
scopic phase coexistence of FM and AF clusters, (iv) For 
some density and A — J5 combination, the ground state 
could be a spin glass. 

Although the phases above can be motivated easily, the 
electrical character of the ground state, or the temper- 
ature dependence of magnetic and transport properties, 
or the response to an applied magnetic field, are still 
not well understood. A comprehensive understanding of 
these effects within the relatively simple model in Eqn.l 
would be the first step in approaching the even richer 
variety of phases in the manganites, where the lattice 
degrees of freedom are also active. This calls for a new 
technique, handling spatial and thermal fluctuations, the 
formation of clusters, and the effect of electron localisa- 
tion. We propose and extensively benchmark such a real 
space technique in this paper. To appreciate the need 
for a new method let us quickly review the current ap- 
proaches to the Hamiltonian above. 

A. Theoretical approaches 

The approaches can be broadly classified into three cat- 
egories. These are: (i) Exact variational calculations 13 at 
T = 0, and generalisation 14-16 toT^O via approximate 
mean field techniques. Let us call these methods varia- 
tional mean field (VMF), for convenience, (ii) Dynamical 
mean field theory (DMFT) based calculations 17 ' 18 which 
map on the lattice model to an effective single site prob- 
lem in a temporally fluctuating medium. Apart from a 
formal limit d — > 00, where d is the number of spatial 
dimensions, there are no further approximations in the 
theory. (Hi) Real space, finite size, Monte Carlo (MC) 
simulations 19-23 of the coupled 'spin-fermion' problem, 
treating the core spin as classical. 

We can set a few indicators in terms of which the 
strength and weakness of various approximations can be 
judged. These are, tentatively: 

f . The ability to access ground state properties. 

2. Ability to handle fluctuations, and accuracy of T c 
estimate. 

3. The ability to access response functions, e.g, trans- 
port and optical properties. 

4. Treatment of disorder effects: Anderson localisa- 
tion and cluster coexistence. 



5. Ability to handle Hubbard interactions, and quan- 
tum effects in spins and phonons. 

6. Computational cost and finite size effects. 

1. Variational calculations 

The variational calculations attempt a minimisation of 
the energy of the (clean) system, at T = 0, with respect 
to a variety of ordered spin configurations. The optimal 
configuration {Si} min for specified J5, [i, etc, is accepted 
as the magnetic ground state. The energy calculations 
are relatively straightforward, since the electron motion 
is in a periodic background. The method has been used 
to map out the ground state phase diagram of DE model 
with AF superexchange in two and three dimension 10 ' 11 . 
The approach, however, can only be approximately im- 
plemented at finite temperature 14-16 . One has to cal- 
culate a spin distribution instead of just targeting the 
ground state, and estimating the energy of an electron 
system in a spin disordered background is non trivial. 
Due to the mean field character of VMF, fluctuation ef- 
fects are lost and transition temperatures are somewhat 
overestimated. The method is focused on thermody- 
namic properties so there is no discussion of transport, 
etc, within this scheme (with one exception 14 ). Disorder 
effects have been included, approximately 14 , in some of 
these calculations. Variational methods can provide in- 
dication of phase coexistence 10 ' 11 at T = 0, or, approxi- 
mately, at finite temperature 15,16 , but cluster coexistence 
in a disordered system is beyond its reach. The method 
has not been generalised to include quantum many body 
effects. Finite size effects in this approach are small and 
the method is relatively easy to implement. 

2. Dynamical mean field theory 

The single site nature of the DMFT approximation be- 
comes exact in the limit of 'high dimensions'. DMFT can 
access both ground state and finite temperature prop- 
erties, but the effective single site approximation can- 
not capture spatial fluctuations, or a non trivial para- 
magnetic phase. The 'mean field' character leads to an 
overestimate of T c , and also the inability to differenti- 
ate between two and three dimensional systems. Being 
a Greens function based theory DMFT can readily ac- 
cess response functions. However, effects like Anderson 
localisation or cluster coexistence, which require spatial 
information, cannot be accessed 24 . The method can han- 
dle many body, Hubbard like, interactions and quantum 
dynamics in all the variables involved, although such cal- 
culations are quite difficult. DMFT is defined directly 
in the thermodynamic limit, so there are no finite size 
effects. The calculations are relatively easy, when quan- 
tum many body effects are not involved, and have been 
a major tool in exploring phenomena in the manganites. 
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The limitations of DMFT become apparent as we con- 
sider the more complicated phases that can arise in our 
model. For instance in the strong disorder problem 7 , 
when there is a possibility of electron localisation, the 
DMFT approach cannot access the insulating phase 24 . 
Neither can it access the spatially inhomogeneous nature 
of freezing, and the persistence of strong spin correlations 
above the bulk T c . Similarly, in the problem of compet- 
ing double exchange and superexchange, in the presence 
of weak disorder, the system breaks up into interspersed 
'ferro-metallic' and 'AF-insulating' regions 8,12 . A com- 
plicated variant of this coexistence effect has been exten- 
sively studied in manganite experiments 25 . The 'single 
site' nature of DMFT cannot access cluster coexistence, 
except possibly in an averaged sense. The transport and 
metal-insulator transitions that can occur in this situa- 
tion also remain inaccessible. So, there are important 
qualitative effects beyond the reach of DMFT, in systems 
where spatial inhomogeneity is important. 



3. Monte Carlo 



suits have been reported 11 ' 23 for the various phases of 
double exchange competing with superexchange antifer- 
romagnetism. No transport results, however, have yet 
been presented within this framework. 

Our method, described in the next section, is devel- 
oped in this spirit. It is a real space Monte Carlo ap- 
proach with the key advantage that it avoids the iter- 
ative TV 3 diagonalisation step. We extract an effective 
Hamiltonian for the core spins from the coupled spin- 
fermion problem, through a self-consistent scheme. We 
can work at arbitary temperature, handle strong dis- 
order, and have better control on 'cluster physics' and 
transport properties due to our significantly larger sys- 
tem size, N ~ 10 3 . 

In the next section we describe our approximation and 
its computational implementation in detail. Following 
that we describe our results on the 'clean' DE model 
in two and three dimension. We will discuss results on 
thermodynamics, spectral features, resistivity and opti- 
cal conductivity, in most of these cases, and compare with 
exact simulation results. We will also highlight systemat- 
ically the size effects in transport and optical properties. 



The finite size real space approach uses the Metropolis 
algorithm to generate equilibrium configurations of the 
spins at a given temperature. Monte Carlo calculations 
on classical systems with short range interactions involve 
a cost 0{zN) for a system update, with z being the co- 
ordination number on the lattice and N the system size. 
In the spin-fermion problem, however, the 'cost' of a spin 
update at a site has to be computed from the fermion free 
energy. If one uses direct diagonalisation of the Hamil- 
tonian to accomplish this, the cost per site is 0(N 3 ), the 
cost for a 'system update' is a prohibitive TV 4 . All this is 
after ignoring quantum many body effects. Current MC 
approaches have not been generalised to handle Hubbard 
like interactions. 

Despite the severe computational cost, this method, 
which we will call ED— MC (exact diagonalisation based 
MC), has been successfully used to clarify several aspects 
of manganite physics, and DE models in general. System 
sizes accessible are ~ 100 at most (recent algorithms 22 
have enhanced this somewhat), with 50 — 60 being more 
typical. This method can provide an outline of the fi- 
nite temperature magnetic phase diagram, reveal major 
spectral features, and even yield the basic signatures of 
cluster coexistence. However, as is obvious from the ac- 
cessible N, the finite size gaps are much too large for any 
reasonable estimate of d.c transport properties, and the 
small linear dimension available, in two or three spatial 
dimension, allows only a preliminary glimpse of coex- 
istence physics. The size limitation apart, the method 
is exact and comprehensive, with none of the problems 
of standard quantum Monte Carlo (QMC). An exten- 
sion of this approach to larger system sizes would al- 
low exploration of several unresolved issues in manganite 
physics. Apart from the ED based MC, 'hybrid MC re- 



II. METHOD 

A. The Ju/t — > oo limit 

We have already written down our basic Hamiltonian 
in Eqn.l. The transformation and projection described 
in the next couple of paragraphs is well known, but we 
repeat them here for completeness. 

Working at large Jujt it is useful to 'diagonalise' the 
J^Si.fTi term first. The electron spin operator is <?i = 
XL/3 c \ a ®ap c iPi wnere the <7^g are the Pauli matrices, 
and this 2x2 problem has eigenvalues ±Jh/2. The 
cigenfunctions are linear combinations of the standard 
'up' and 'down' z quantised fermion states at the site: 
l\n = E Q ^a c L- Thc lower energy state, 7^, a linear 
combination of thc form A^c^ + A\ 2 c\^ is at energy 
—Jh/2 and has fermion spin parallel to the core spin Sj. 
The orthogonal linear combination, 7^, has fermion spin 
anti-parallel to the core spin and is at energy +J#/2. 
The amplitudes A 1 are standard 9 . 

In thc 7 basis, the Hunds coupling term becomes 
-(Jj?/2)(7| ; 7 ii - 7] u 7i„) at all sites. The intersite 
hopping term, however, picks up a non trivial depen- 
dence on nearest neighbour spin orientation, tijc\ a Ci a — > 

E Q /3 tijdiflL'yjP where a, f3 refer to thc it, I indices. g"f 
arises from the product of the two transformations at site 
i and site j, and we will describe its specific form later. 
Since the canonical transformation is local, the density 
operator J2 a c\ a c ia -> (7^ + 7i„7m)- 

At finite Ju/t this is just a transformation from the 
'lab frame' to a local axis and the 'up' and 'down' spin 
fermions get mapped to (l,u), but we still have to solve 
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a mixed 'two orbital problem'. However, if Ju/t — > oo 
then all the 'anti-parallel' 7/jO) states get projected out 
and we can work solely in the subspace of states created 
by 7^ . In this space, the Hamiltonian assumes a simpler 
form: 

H e i = -t^2( 9iilhj + h.c )+ ^(e; - ii)n t 

(ij) ' 

= -t E fa ( ^ ihi + ^) + E( e «- ^ ( 2 ) 

(ij) 1 

where we have dropped the superfluous 11 label in gij, 
and absorbed — J#/2 in the chemical potential. The 
hopping amplitude = fije % ^ ij between locally aligned 
states, can be written in terms of the polar angle 

ft 9 

and azimuthal angle (</>i) of the spin Si as, cos^-cos-^- 
+sin^sin^e~ l It is easily checked that the 

'magnitude' of the overlap is = + Sj.Sj)/2, while 
the phase is specified by tan&ij — Im(gij)/ Re(gij). 

This problem can be viewed as a quadratic 'spin- 
less fermion' problem with core spin dependent hopping 
amplitudes. The fermions move in the background of 
quenched disorder e, and 'annealed disorder' in the {Si}, 
where the second brackets indicate the full spin configu- 
ration. To exploit the nominally 'non interacting' struc- 
ture of the fermion part we need to know the relevant 
spin configurations, {Si}, or, more generally, the distri- 
bution P{Si}, controlling the probability of occurence of 
a spin configuration. 

B. Effective Hamiltonian for spins 

The partition function of the system is Z = 
JVSiTre-? 11 . To extract P{Si} note that for a sys- 
tem with only spin degrees of freedom, Z will have the 
form J VSie~P H { s ^ . Comparing this with the partition 
function of the spin-fermion problem we can use 

J VS t Tre-P H = J VS^-^f^ 

from which it follows that 

H eff {S i } = ~logTre-' 3H 

PiS,} oc e -^/H s *> (3) 

The trace is over the fermion degrees of freedom. In our 

case 

H ef/ = -hogTre-^ + J s ]T S i .S j (4) 

P (ij) 

The principal difficulty in a simulation, and quite gen- 
erally in spin-fermion problems, is in evaluating the first 
term on the r.h.s above for an arbitrary spin configura- 
tion. This is the origin of the iV 3 factor in the exact MC. 



Our key proposal, whose analytic and numerical justifi- 
cation we provide later, is 

-hogTre-^n-^Diih (5) 

P (ij) 

where is an effective 'exchange constant' to be deter- 
mined as follows. Define the operator = (e 1 * 4 -^^ + 
h.c). This enters the 'hopping' part of the electron 
Hamiltonian. In any specified spin configuration {/, $} 
we can calculate the correlation function Dij{f, $} = 
Z~ t TrTije~P Hel , where Z e i is the electronic partition 
function in the specified background. The exchange 
that finally enters H e fj is the average of Dij{f, $} 
over the assumed equilibrium distribution, i.e: = 
J VfT>$P{f, $}Dy {/, $} where we denote a spin config- 
uration interchangeably by {/, $} or {S}. Qualitatively, 
the 'effective exchange' is determined as the thermal av- 
erage of a fermion correlator over the assumed equilib- 
rium distribution. Let us bring together the equations 
for ready reference. 

H e i = -t E fij^ij + /~2( e i - 

{ij) 1 

hi = ^(l + S,.S,)/2 
H eff {S} = -hogTre-^ H " + J S ^2 S ^ S J 

1 m 

~ - E D ijfv + J S E Si ' S j 

(ij) (ij) 

Dij = (<Ph))h.„ (6) 

The ED— MC approach 'solves' for physical properties 
by using the first four equations above: equilibriating 
the spin system by using H e ff, which itself involves a 
solution of the Schrodinger equation for the electrons. 

Our method approximates the 'exact' H e ff by the form 
specified in the fifth equation and computes an exchange, 
rather than equilibrium configurations themselves, by 
fermion diagonalisation. The sixth equation indicates 
how the 'loop' is closed. We will refer to this method as 
"Self Consistent Renormalisation" (SCR) 26 , or the H e ff 
scheme. 

The nonlinear integral equation for the is solved to 
construct the 'classical spin model' for a set of electronic 
parameters, disorder realisation, and temperature. Al- 
though the assumption about H e ff seems 'obvious', and 
in fact something similar, but simpler, had been explored 
early on by Kubo and Ohata 4 , and recently by Calderon 
and Brcy 20 , the power of the method becomes apparent 
in disordered systems or in the presence of competing 
interactions. In these cases the solutions Dy can be spa- 
tially strongly inhomogeneous, and dramatically temper- 
ature dependent. The properties of such systems are far 
from obvious. 
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The equilibrium thermal average of any fermion oper- 
ator, or correlation function, O, can now be computed 
using the self-consistent distribution as: 

((6)) = J V{S}P{S}0(S) (7) 

The average O(S) is computed on a spin configuration 
{S,}, with the configurations themselves picked accord- 
ing to the effective Boltzmann weight oc e~P H *ss . 

We have not written the equation for fi. Since we would 
typically want to work at fixed density rather than fixed 
chemical potential, we employ the procedure above to 
calculate n and iterate \x till the 'target' density is ob- 
tained. In actual implementation, discussed later, the \x 
'loop' and the 'loop' run simultaneously. We next 
discuss the analytic underpinning of our method before 
moving to numerical results. 



C. Analytic limits 

The central problem in DE models is construction of an 
energy functional for arbitrary spin configurations {/,$}. 
This information is contained in the fermion free energy, 
—TlogTre^ 1311 ' 1 as we have seen. We study two limits 
below, where the leading effects are well captured by our 
effective Hamiltonian. 



1. Low temperature 

If we ignore disorder and AF coupling, for simplicity, 
and if the free energy of the fermions can be approxi- 
mated by the internal energy, then Dij{f, $} contains the 
necessary information about the energy of any spin con- 
figuration: £{/,<&} = H eff {f,*} V !),,{/. <!-[/,,. 
The configuration dependent correlation function, how- 
ever, is hard to calculate, since it requires a solution of 
the Schrodingcr equation for each spin configuration. 

At low temperature, as the spins gradually ran- 
domise, the system explores configurations {/, $} near 
the ground state in the energy landscape. The relevant 
Aj{/, $} ~ D^+SDij{f, $}, where D° is the 'exchange' 
computed on the ground state, and SD is the variation. 
At low T, such that the relevant SD <C D°, we can ne- 
glect the variation, SD, between configurations, and the 
'effective Hamiltonian' assumes the form: 

lim T ^ H eff ~ - £ V%h = - E D % V( l + Si.SjO/2 

(ij) (ij) 

As we will see in the simulations this approximation is 
remarkably good in the simple DE model almost upto 
T c /2. At higher T the 'renormalisation' of D becomes 
important. 



2. High temperature 

For T c /T <C 1, cumulant expansion yields an asymp- 
totically exact effective Hamiltonian: 

H eff ~ limpt^a- ^lnTr(l+(3H + ^p 2 H 2 + ..) 

The leading contribution from this is: 

Hffi T ~-n(l-n)/3i 2 £/ 2 . 

(ij) 

This apparently has a structure different from that of our 
H e ff, and additionally an 'effective coupling' falling off 
as 1/T. In fact our coupling D has the same form, as can 
be checked by evaluating ((Tij)) in a high temperature 
expansion. This quantity also depends on n(l — n), to 
allow hopping, and falls off as 1/T since it is non local. 
The self- consistent calculation of the effective exchange, 
now based on the high temperature phase rather than 
the ground state, ensures that the leading contribution 
to the energy is well captured. The physical consequence 
of the 1 /T effective exchange is that the susceptibility of 
the DE model does not have the Curie- Weiss form that 
one expects for Heisenberg like models 27 . 

The next order in series expansion will generate terms 
of the form: 

J2f ij f j kfkifue i ^ + ^ + -\ 

ijkl 

summed over the minimal plaquette. Higher powers in 
jit involve longer range excursion of the fermions, but 
the limited data available from exact simulations suggests 
that the critical properties of double exchange are similar 
to that of short range spin models. 

Although the procedure above can be extended to ex- 
tract an 'exact' effective Hamiltonian to high order in 
fit, we know of no such attempt. The only series expan- 
sion results available are on the S = 1/2 model, directly 
calculating thermodynamic properties 28 . 

D. Monte Carlo implementation 

Since the ground state of the system is often not known 
it is usual to start from high temperature and follow the 
sequence below in generating the effective Hamiltonian 
and studying equilibrium properties. 

(i) We start at high temperature, T » T c , assuming 
some Dfj (T), where n is the iteration index, and 'equilib- 
riate' the system with this assumed effective Hamiltonian 
(not yet self-consistent), (ii) We compute the average 
((e'* y 7j7j + h.c)) over these (pseudo) equilibrium con- 
figurations. This generates D™ +1 (T). {Hi) Compare the 
generated exchange with the assumed exchange at each 
bond. Accept if within tolerance. If converged, then Dij 
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represents the correct 'exchange' at that temperature. 
Else, replace £>™ by D"^ 1 - (iv) At each temperature 
and iteration, adjust \x as necessary to keep n constant. 

At convergence fermion properties can be calculated 
and averaged over equilibrium MC configurations of the 
spin model. For a disordered system (A ^ 0), the ther- 
mal cycle above has to be repeated for each realisation 
of disorder. In the clean problem, translation invariancc 
forces the exchange to be uniform at all bonds, while 
A ^ generates a bond disordered spin model. 

The computational effort needed in the ED— MC ap- 
proach is oc Nmc x iV 4 , at each temperature, where Nmc 
is the number of MC sweeps (10 3 — 10 4 ), and N the size 
of the system (actually the Hilbcrt space dimension). 
As we have mentioned before, current resources allow 
N m ax ~ 100. Within our H e ff scheme the MC config- 
urations are generated using a short range spin model, 
with cost 0(N). The actual cost is in determining the 
exchange: this is cx Nu er x N av x N 3 , where Ni ter is the 
number of iterations needed to get a converged solution, 
with ~ 10% accuracy per bond, and N av is the averaging 
needed per iteration for generating a reasonable 'equilib- 
rium average'. Typically Ni ter ~ 4 and N av ~ 50. 

We can roughly compare the computational cost of 
ED-MC with the H eff scheme. For ED-MC, the time 
required is, tjv ~ Nmc x TV 4 at a given temperature. 
For the H e ff scheme, tn ~ Nu er x N av x N 3 . Putting in 
the numbers, if resources allow N ~ 100 for the ED— MC 
approach, the same resource will allow N ~ 1000 within 
the H e ff scheme. In terms of computation time, H e ff 
is no more expensive than standard 'disorder average' in 
electronic systems. 



E. Physical properties at equilibrium 

The major physical properties we compute at equilib- 
rium are optical conductivity and d.c resistivity, the den- 
sity of states (DOS), and the magnetic structure factor. 

(i) We estimate the d.c conductivity, cr^c, by using the 
Kubo-Greenwood expression 29 for the optical conductiv- 
ity. In a disordered non interacting system we have: 

The constant A = (ne 2 ) /Tia®. The matrix element 
lap = (ipabxlipp) an d we use the current operator 

jx = ia-oY<i, 7 { c \+xa ,a c i,a ~ h - c )- Tnc etc are sin- 
gle particle eigenstates, for a given equilibrium configu- 
ration, and e a , ep are the corresponding eigenvalues. The 
n a = 8(p — e a ), etc, are occupation factors. 

The conductivity above is prior to thermal or disor- 
der averaging. Our simulations are in a square or cube 
geometry with periodic boundary condition. Given the 
finite size, the S function constraint in o~(u>) cannot be 
satisfied for arbitrary w. We use the following strategy: 



(i) calculate o~i n t(uj) = Jq o~(u)')du)' , at three equispaced 
low frequency points, u>i,u>2,^>3, by summing over the 
delta functions in the appropriate range., (ii) thermally 
average the o~i n t(uj) over the equilibrium configurations, 
(Hi) invert: calculate a numerical derivative via three 
point interpolation, implementing a(u>) = dai n t(uj) / duj. 
The 'bar' on a indicates thermal average. What we call 
the 'd.c. resistivity' is actually the inverse of a low fre- 
quency optical conductivity, computed by the method 
above. We systematically check the stability of our re- 
sults by repeating the calculation for a sequence of system 
size (and reducing oj\, u>2, ^3 accordingly). For N ~ 1000, 
the 'd.c' conductivity is actually computed at u> ~ 0.06. 

Our transport calculation method and some bench- 
marks will be discussed in detail elsewhere 30 . To con- 
vert to 'real' units, note that our conductivity results are 
in units of (ire 2 )/hao. Since the Mott 'minimum' metal- 
lic conductivity, in three dimension, is ~ (0.03e 2 /?iao), 
a = 1 on our scale roughly corresponds to lQ 2 o~Mott- 
The full <j(lo) is computed by computing a int (uj) defined 
above, thermal average, and inversion. 

(ii) Each equilibrium magnetic configuration leads to 
a 'DOS' of the form ^ Q S(u) — e a ), where e a are the single 
particle eigenvalues in that background. The thermally 
averaged DOS that we show involves a Lorentzian broad- 
ening of each S function, as indicated below. 

The sum runs over the eigenvalues obtained in any spin 
configuration, and summed over equilibrium configura- 
tions. We use r ~ 0.1 in our results, although much 
smaller T would still give a smooth spectra at high T. 
(Hi) The magnetic structure factor is calculated as 

€q eq ij 

where i, j run over the entire lattice, and the outer aver- 
age is over equilibrium configurations. 



III. RESULTS 

In this section we provide a comprehensive compar- 
ison of results based on the 'exact' scheme (ED— MC) 
and our effective Hamiltonian approach, for the 'clean' 
DE model, and and extend the study to large sizes using 
the H e ff scheme. Most of our results are on three dimen- 
sional systems, where the simulations are more difficult 
and the results physically more relevant, and we show 
only limited data in two dimensions. The model is trans- 
lation invariant, there are no competing interactions, and 
the low temperature phase is a ferromagnet. 
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A. Magnetism and thermodynamics 

We begin with a comparison of the magnetisation, 
m(T), obtained via ED-MC and SCR on 8 x 8 lattices in 
2d, and 4 3 systems in 3d, with periodic boundary condi- 
tion in all directions. Fig.l compares the m(T) obtained 
via the two schemes at three electron densities. 

Note at the outset that both the DE model and our 
H e ff are 0(3) symmetric and are not expected to have 
long range order at finite T in 2d (in an infinite system) . 
However, as has been demonstrated in the case of the two 
dimensional classical Heisenberg model 31 , 0(3) models 
have exponentially large correlation length at low tem- 
perature in 2d. For a nearest neighbour classical Heisen- 
berg model with |S,| = 1, and exchange J, the low T 
correlation length £(T) - 0.02e 27rJ / T . So, for T < J 
even large finite lattices would look 'fully polarised' and 
one would need to access exponentially large sizes to see 
the destruction of long range order. 

This allows us to define a (weakly size dependent) 
'characteristic temperature' T c h{n) for the 2d DE model 
which marks the crossover from paramagnetic to a nom- 
inally 'ordered' phase. The true ordering temperature of 
strongly anisotropic DE systems, e.g, the layered man- 
ganites, which the planar model mimics, would be de- 
termined by the interplane coupling, but the in plane 
transport would be dictated mainly by the 2d fluctua- 
tions, as here. 

The difference between ED-MC and SCR results in 
2d, Fig.l. (a), is most prominent at the highest density, 
n = 0.41, where the T c h inferred from these small size 
calculations differ by 15% — 20%. At lower density the 
difference is still visible but much smaller. We have indi- 
cated the T c h scales inferred from the two schemes in the 
inset in Fig.l. (a) The difference between the two schemes 
is usually largest in clean high density systems, as we will 
see also in the three dimensional case. However, over the 
entire density range, the maximum deviation is ~ 20%. 
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bols) -vs- SCR results (dotted lines). 
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(b). internal energy, (c). chemical potential and (d). band 
edge. Displayed value is actual value + shift. System size 
4x6x4. Open symbols: ED-MC, dotted lines: SCR. 

Notice that at all n, the low temperature m(T) ob- 
tained via H e ff corresponds almost exactly with results 
based on ED-MC. This works upto - T ch /2. The high 
temperature result within the two schemes is also in close 
correpondence but that is better illustrated in the ther- 
modynamic data, Fig. 2, which we will discuss later. 

Fig.l.(b) shows the results on magnetisation in the 
three dimensional problem at three densities, comparing 
results based on ED— MC and H e ff. As in two dimen- 
sion the difference in the estimated T c is greatest near 
the band center, being ~ 15% — 20%, the correspondence 
improving as we move to n < 0.2. As before, the exact 
and approximate m(T) match at low T for all densities. 

Fig. 2 which shows the thermodynamic indicators in the 
3d case reveals that itself is virtually indistinguish- 
able in the two schemes. The correlation = ((1%-)} 
can be evaluated as an equilibrium average in an exact 
simulation also, although there it does not feed back into 
the calculation. The match between the Z)'s computed 
in two different schemes, and across the density range, 
suggests that the difference in m(T) seen near half-filling 
is not due to different numerical values of D, but the 
assumed form of H e ff. We either need a more sophis- 
ticated definition of the finite temperature D, or a dif- 
ferent form of H e f f to bring the high density results of 
H e ff in closer correspondence with ED— MC. Notice that 
the .D's are only weakly temperature dependent and the 
m(T) at low temperature could have been obtained by 
setting D(T) — D(Q). In fact over the temperature range 
— T c the qualitative physics can be accessed without the 
thermal 'renormalisation' of the exchange. However, for 
T » T c the renormalisation is important, as suggested 
earlier by the high temperature expansion. 

The results on all thermodynamic indicators, D(T), 
E(T), n(T) and E b (T), Fig.2, show the close corre- 
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spondence between results of the exact and approximate 
scheme. The -D's are almost temperature independent 
in the range — T c and hardly distinguishable between 
ED— MC and H e ff, suggesting that effects beyond our 
effective Hamiltonian —DJ2fij is needed to accurately 
describe the magnetic transition at the band center. The 
overall behaviour is similar in 2d as well so we are not 
presenting the 2d data. 

We extend the H e ff scheme to large system size, and 
study the magnetism in 32 2 and f 3 lattices. Fig. 3 shows 
the results on m(T), and the inset shows the T c inferred 
from these simulations. The maximum T c , occuring at 
band center is ~ 0.2t which, with t ~ (100 — 150)mcV, 
will be in the range 200— 300K. These numbers are typical 
of high electron density Hunds coupled systems, and are 
in the right ballpark when compared to the manganites 32 . 
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FIG. 3. Magnetisation based on H e f f in (a). 2d with 30x30 
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acteristic temperature scales inferred from m(T) 
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FIG. 4. Thermodynamic properties in the 3d case com- 
puted with H e ff, system size N = 10 x 10 x 10. 



Fig. 4 shows the thermodynamic indicators computed 
within the H e ff scheme on 10 3 in 3d. The strong temper- 
ature dependence in fj, and E^, seen also at small sizes, 
arise from the 'band narrowing' effect of spin disorder 
which reduces the mean hopping amplitude with increas- 
ing temperature. 



B. Density of states 

Fig. 5 shows the density of states (DOS) computed at 
n = 0.3, four temperatures, and for a small, 4 3 , and 
a large, 10 3 , system. The mean level spacing at high 
temperature (where the spins are completely disordered) 
is ~ 12/i 3 which is ~ 0.01 at L = 10 and ~ 0.18 at L = 
4. For T — > 0, the polarised ferromagnetic state leads 
to large degeneracy and the level spacings could be more 
than 10 times larger than the high temperature value. 
We have broadened all 5 functions by T — 0.1, so that 
the high temperature L = 4 spectra looks reasonable. 
With this broadening the L = 10 data looks reasonable 
even below T c . 

This comparison highlights the unreliability of small 
size data in inferring spectral features over most of the 
interesting temperature range. Small sizes can often pro- 
vide reasonable results on energetics, but on spectral fea- 
tures and, more importantly, on low frequency transport, 
they are completely unreliable. 



C. Optical properties 

Fig. 6 shows the optical conductivity, a(cj). The opti- 
cal conductivity is a vital probe of charge dynamics in 
the system. Our data in the main panel, Fig. 6, is for a 
8x8x8 geometry. At the lowest temperature there is an 
artificial 'hump' in <j(ui) which we think arises because 
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the polarised three dimensional system has large degen- 
eracy, and finite size effects are stronger than in two di- 
mension. Nevertheless, there are some notable features in 
a(cu), (i) the conductivity is Drude like, (ii) there is rapid 
reduction in low frequency spectral weight with increas- 
ing temperature, with some transfer to high frequency 
(Hi) the weight in a(u>) is not conserved with increas- 
ing temperature, the loss is related to the suppression of 
kinetic energy with increasing spin disorder. 




FIG. 6. Optical conductivity based on H e ff in three di- 
mension. System size 4x4x4 (inset) and 8x8x8 (main 
panel), density n — 0.3. Symbols in the inset are same as in 
the main panel. 
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FIG. 7. (a) — (b). Magnetisation and normalised resistiv- 
ity at n = 0.3 in (a). 2d and (b). 3d. Insets to (a), and 

(b) . show the size dependence in the resistivity (see text). 

(c) . — (d). Density dependence of p(T), in 2d and 3d respec- 
tively, with system sizes 32 x 32 and (d) 10 x 10 x 10. 
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D. Resistivity 

Finally, we look at the resistivity, which, surprisingly, 
has seen little discussion. Fig. 7 shows the correlation be- 
tween the fcrromagnet to paramagnet transition and the 
rise in p(T). We have normalised p(T) in Fig.7.(a)-(b) by 
the value at T = 0.4. The 'absolute' resistivity is shown 
in Fig.7.(c)-(d). Unlike mean field treatments which treat 
the paramagnetic phase as completely 'uncorrected' and 
would yield a 'flat' resistivity for T > T c h (or T c ), there is 
a significant increase in p(T) with rising temperature in 
the 'paramagnetic' phase as the short range spin correla- 
tion is gradually lost and the system heads towards the 
fully spin disordered phase. The general rise in p(T) in 
the paramagnetic phase happens in both 2d and 3d, but 
surprisingly in 2d most of the rise seems to occur after 
the drop in m(T), rather than across T c as one sees in 
three dimension. 

For a check on the reliability of the computed p(T) the 
inset in Fig. 7. (a) shows the 'resistivity' computed on L x 
L geometry for L = 8, 16, 32 across the full temperature 
range. The L = 8 result has the same problem that we 
discussed in the context of cr(oj). The system essentially 
behaves as an 'insulator' at low T due to the finite size 
gap. The L — 16 data has similar upturn, but at a 
lower temperature. The data at L = 24 (not shown) and 
L = 32 are stable down to T ~ 0.02 and almost coincide, 
suggesting that except at very low temperature, results 
on these sizes are representative of bulk transport. 

The resistivity in the 3d case differs from 2d in that 
the major rise in p(T) occurs around T c in the 3d case, 
while it occurs beyond T c h in the 2d case. Fig.7.(b) shows 
m(T) correlated with the normalised p(T), and the rise 
is reminiscent of the Fisher-Langer result 33 in weak cou- 
pling electron-spin systems. The inset in Fig.7.(b) shows 
the stability of the transport result in 3d for L > 8, and 
the unreliability for L ~ 4. 

Fig.7.(c)-(d), shows the absolute resistivity for a few 
densities. The 'high temperature' 3d resistivity at T ~ 
3T C is approximately 15 — 25, in the density range shown, 
which in real units would be ~ (1 — 2)mOcm, roughly the 
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high T resistivity of Lai-^Sr^MnOs for x > 0.4. 

Fig. 8. shows the correlation between dp/dT and the 
specific heat in 2d and 3d. Above T c and in 3d, panel 
(b), dp/dT seems to match Cy very well, as expected 
from the perturbative results of Fisher and Langer 33 . In 
2d however the correspondence is poor, probably due to 
incipient localisation effects in the resistivity. For T < T c , 
even in 3d, the behaviours of Cy and dp/dT are different 
because the rise in m(T) affects the scattering rate, as is 
already known 33 . 

The validity of the 'weak coupling' results of Fisher- 
Langer, originally illustrated for a Heisenberg model, in 
this 'strong coupling' spin-fermion system may seem sur- 
prising. There are two reasons why the correspondence 
holds here: (i) the resistivity in the DE model arises from 
spin disorder induced weak fluctuations in the hopping 
amplitude, and is in the perturbative regime, and (ii) our 
magnetic model, H e f f , is effectively short range, and the 
critical properties of spin fluctuations are the same as in 
the Heisenberg model. 



IV. CONCLUSION 



In this paper we proposed a new Monte Carlo tech- 
nique that allows access to large system sizes but retains 
the correlated nature of spin fluctuations in the double 
exchange model. Combining this MC technique with a 
transport calculation based on the exact Kubo formula 
we presented a comprehensive solution of the model, in- 
cluding magnetism, thermodynamics, spectral features, 
transport, and optics. 

This paper benchmarked the scheme for the clean dou- 
ble exchange model, where the complicated consistency 
and thermal renormalisation involved in the scheme are 
not crucial for a qualitative understanding. However, 
when we move to disordered systems 7 , or non ferro- 
magnetic ground states, or the regime of multiphase 
coexistence 8 , the full power of a 'bond disordered' ef- 
fective Hamiltonian, with non trivial spatial correlation 
between the bonds, becomes apparent. 

For the clean ferromagnetic case one may try to 
improve the self-consistency scheme to obtain better 
correspondence 20 ' 34 with ED— MC results. However, 
given the complexity of the current scheme, and the range 
of possibilities that it offers, we think it is more important 
to exploit the present scheme to resolve the outstand- 
ing qualitative issues first. Finally, although the entire 
scheme is presently implemented numerically, it would 
be useful to make analytic approximations within this 
framework to create greater qualitative understanding. 

We acknowledge use of the Beowulf cluster at HRI. 
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